#############################################
# This file contains user-written functions #
# for use in replication of:                #
#                                           #
# "The Power to Hurt and the Effectiveness  #
# of International Sanctions"               #
#                                           #
# By Kerim Can Kavakli,                     #
# J. Tyson Chatagnier,                      #
# and Emre Hatipoglu                        #
#                                           #
# Forthcoming in the Journal of Politics    #
#                                           #
# File created 25 January 2019              #
#############################################

####################
# Define functions #
####################

localLogitLik <- function(
                          params, 
                          X, 
                          Y, 
                          q1, 
                          q2, 
                          q3, 
                          nclus        = 4,
                          cluster.init = FALSE,
                          quietly      = FALSE
                          ) {
  require(matrixStats)
  if (cluster.init == FALSE) {
    require(doParallel)
    require(foreach)
    cl <- makeCluster(nclus)
    registerDoParallel(cl)
  }
                        # Telling R to use a cluster of size nclus.
                        # If this is inside something else (as is
                        # the case with cross-validation), then 
                        # we can set the flag to TRUE, in order 
                        # to keep from doing this many times.
  if (any(is.na(params)) | length(params) < 3) {
    stop("'params' must be a vector of length 3.")
  }
                        # Smoothing parameter values must be
                        # entered as vectors of length 3. Missing
                        # values are not allowed.
  h <- params[1]
  delta <- params[2]
  lambda <- params[3]
  if (length(unique(X[, 1])) == 1) {
    X.const <- X
    X <- X[, 2 : ncol(X)]
  } else {
    X.const <- cbind(
                     1,
                     X
                     )
  }
  mat.b <- foreach(
                   i         = 1 : nrow(X), 
                   .combine  = "rbind", 
                   .packages = "matrixStats", 
                   .export   = c(
                                 "contWt", 
                                 "kernelWeight", 
                                 "ordWt", 
                                 "unordWt"
                                 )
                   ) %dopar% {
                     kernel.wt <- (
                                   contWt(
                                          q1, 
                                          X, 
                                          i, 
                                          h
                                          ) * 
                                   ordWt(
                                         q1, 
                                         q2, 
                                         X, 
                                         i, 
                                         delta
                                         ) * 
                                   unordWt(
                                           q3, 
                                           X, 
                                           i, 
                                           lambda
                                           )
                                   )[-i]
                        # Calculate the kernel weight for each observation.
                        # foreach allows parallel processing.
                       if (mean(kernel.wt) > 0) {
                        # This if statement is meant to take care of the
                        # situations in which we get a kernel weight of
                        # zero on all observations, leading to the
                        # 'fit not found' error.
                         glm.fit(
                                 X.const[-i,], 
                                 Y[-i], 
                                 weights = kernel.wt, 
                                 family  = binomial(link = "logit")
                                 )$coef
                       } else {
                         rep(
                             NA, 
                             times = ncol(X.const)
                             ) 
                       }
    }
                        # Estimate leave-one-out fitted likelihood for 
                        # each observation, using the kernel weight
                        # obtained above.
                        # Use the coefficients to calculate likelihoods
                        # below.
  if (cluster.init == FALSE) {
    stopCluster(cl)
  }
                      # Turn off the cluster if we initialized it here.
  xb <- rowSums(X.const * mat.b)
                        # This performs the same operation as 
                        # applying '%*%' to each row.
  lnp0 <- log(1 / (1 + exp(xb)))
  lnp1 <- log(1 / (1 + exp(-xb)))
                        # Calculate probabilities of 0 and 1, 
                        # respectively.
  lnp0 <- ifelse(
                 lnp0 == -Inf, 
                 -xb, 
                 lnp0
                 )
  lnp1 <- ifelse(
                 lnp1 == -Inf, 
                 xb, 
                 lnp1
                 )
                        # These two lines eliminate the problem of 
                        # machine zero and rounding error.
                        # The log of (1 + exp(x)) approaches x 
                        # as x goes to infinity, so at the values
                        # where this is a problem, the change will 
                        # be negligible, and should be
                        # superior to what R is trying to do.
  lik <- Y * lnp1 + (1 - Y) * lnp0
  log.lik <- sum(
                 lik, 
                 na.rm = TRUE
                 )
                        # Calculate the log-likelihood.
  log.lik <- ifelse(
                    log.lik == 0, 
                    -1e20, 
                    log.lik
                    )
                        # If we run into problems of all NAs,
                        # set the log-likelihood to a very 
                        # large negative value. This ensures
                        # that we don't accidentally pick too
                        # small smoothing parameters.
  if (!quietly) {
    cat(
        "\n", 
        "Parameters: ", 
        "h =", 
        h, 
        ", delta =", 
        delta, 
        ", lambda =", 
        lambda, 
        "\n", 
        "Log-Likelihood: ", 
        log.lik, 
        "\n"
        )
                        # Print the results of the estimation,
                        # if desired.
    }
    return(log.lik)
}

crossValidation <- function(
                            stval.mat, 
                            X, 
                            Y, 
                            q1, 
                            q2, 
                            q3,
                            nclus     = 4,
                            bounds.lo = c(
                                          3, 
                                          0.01, 
                                          0.01
                                          ) 
                            ) {
                        # Note that stval.mat must be a 
                        # three column matrix, but can have
                        # any (positive) number of rows.
  require(doParallel)
  require(foreach)
  data.set <- cbind(
                    Y,
                    X
                    )
  data.nona <- na.omit(data.set)
  Y <- data.nona[, 1]
  X <- data.nona[, -1]
                        # Can't handle NAs, so eliminate those.
  params.mat <- matrix(
                       nrow = nrow(stval.mat), 
                       ncol = ncol(stval.mat)
                       )
  llik.vec <- vector(length = nrow(stval.mat))
  cl <- makeCluster(nclus)
  registerDoParallel(cl)
  for (i in 1 : nrow(stval.mat)) {
      stval.bw <- stval.mat[i,]
                        # Using start values from the matrix
      bw.est <- optim(
                      par     = stval.bw, 
                      fn      = localLogitLik,
                      hessian = F, 
                      method  = "L-BFGS-B", 
                      lower   = bounds.lo, 
                      upper   = c(
                                  Inf, 
                                  1, 
                                  1
                                  ),
                      control = list(
                                     fnscale = -1, 
                                     trace   = 1, 
                                     REPORT  = 1, 
                                     maxit   = 2000, 
                                     factr   = 1e-300
                                     ), 
                      X  = X, 
                      Y  = Y,
                      q1 = q1, 
                      q2 = q2, 
                      q3 = q3,
                      nclus = nclus,
                      cluster.init = TRUE
                      )
                      # The coefficients returned here should be the 
                      # optimal bandwidth parameters.
                      # Need to put upper and lower bounds on the
                      # parameters, as they are bounded by definition.
                      # Can adjust bounds.lo to small positive numbers in
                      # order to avoid problems of near-zero weighting.
                      # Setting cluster.init to TRUE because we are
                      # initializing within the function.
      params.mat[i, ] <- bw.est$par
      llik.vec[i] <- bw.est$value
  }
  stopCluster(cl)
  return(
         list(
              par = params.mat, 
              llik = llik.vec
              )
         )
}

localLogit <- function(
                       X.mat, 
                       Y.vec, 
                       bw.par, 
                       profile.mat,
                       nclus = 4, 
                       q1, 
                       q2, 
                       q3, 
                       tval  = 1.96
                       ) {

  require(matrixStats)
  require(doParallel)
  require(foreach)
  data.set <- na.omit(
                      cbind(
                            Y.vec, 
                            X.mat
                            )
                      )
                        # Remove missing data.

  Y <- data.set[, 1]
  X <- data.set[, -1]
  k <- ncol(X)
  n <- nrow(profile.mat)
  h <- bw.par[1]
  delta <- bw.par[2]
  lambda <- bw.par[3]
  kernel.wt <- matrix(
                      nrow = n,
                      ncol = length(Y)
                      )
  for (i in 1 : n) {
    kernel.wt[i, ] <- contWt(
                             q1, 
                             X, 
                             i, 
                             h, 
                             profile.mat
                             ) * 
                      ordWt(
                            q1, 
                            q2, 
                            X, 
                            i, 
                            delta, 
                            profile.mat
                            ) * 
                      unordWt(
                              q3, 
                              X, 
                              i, 
                              lambda, 
                              profile.mat
                              )
  }
  cl <- makeCluster(nclus)
  registerDoParallel(cl)
  list.mod <- vector(
                     "list",
                     length = n
                     )
  list.mod <- foreach(i = 1 : n) %dopar% {
                glm(Y ~ X, 
                    weights = kernel.wt[i, ], 
                    family  = binomial(link = "logit")
                    )
  }
  stopCluster(cl)
                        # This produces a list of estimated models
                        # with n different entries. It should be
                        # a separate model entry for each
                        # profile.
  betas <- t(
             sapply(
                    list.mod, 
                    coef
                    )
             )
                        # This extracts the coefficients from each of
                        # the models within list.mod, and turns them
                        # into an n x k matrix.
  mat.c <- cbind(
                 1, 
                 profile.mat
                 )
  xb <- rowSums(
                mat.c * betas, 
                na.rm = T
                )
  pred.probs <- plogis(xb)
                        # This calculates predicted probabilities for
                        # each profile of values (row of profile.mat).
  term.1 <- (pred.probs ^ 2) * (1 - pred.probs) ^ 2
  se.vec <- vector(length = n)
  for (i in 1 : n) {
    vhat <- calcVar(
                    X    = X,
                    Y    = Y,
                    beta = list.mod[[i]]$coef,
                    wt   = kernel.wt[i, ]
                    )
    term.2 <- t(as.matrix(mat.c [i, ])) %*% vhat %*% as.matrix(mat.c[i,])
    se.vec[i] <- sqrt(term.1[i] * term.2)
  }
                        # This calculates standard errors for the
                        # predicted probabilities.
  ci.lo <- plogis(xb) - (tval * se.vec)
  ci.lo <- ifelse(
                  ci.lo < 0,
                  0,
                  ifelse(
                         ci.lo > 1,
                         1,
                         ci.lo
                         )
                  )
  ci.hi <- plogis(xb) + (tval * se.vec)
  ci.hi <- ifelse(
                  ci.hi < 0,
                  0,
                  ifelse(
                         ci.hi > 1,
                         1,
                         ci.hi
                         )
                  )
                        # This gives confidence intervals based on
                        # the predicted probabilities and standard errors.
  return(
         list(
              b       = betas,
              p       = pred.probs, 
              profile = profile.mat, 
              lo      = ci.lo, 
              hi      = ci.hi
              )
         )
}

calcVar <- function(
                    X,
                    Y,
                    beta,
                    wt
                    ) {
  if (any(X[, 1] != 1)) {
    X <- cbind(
               1, 
               X
               )
  }
  Xb <- X %*% beta
                        # This gives us a vector of
                        # Xb values for each point in
                        # the dataset, given a particular
                        # set of betas.
  lambda <- plogis(Xb)
  lambda.bar <- 1 - lambda
  term.outer <- 0
  term.inner <- 0
  for (i in 1 : nrow(X)) {
    term.outer <- term.outer + (lambda[i] * lambda.bar[i] * X[i, ] %*% t(X[i, ]) * wt[i])
    term.inner <- term.inner + ((Y[i] - lambda[i]) ^ 2 * X[i, ] %*% t(X[i, ]) * wt[i] ^ 2)
  }
  var.x <- solve(term.outer) %*% term.inner %*% solve(term.outer)
  return(var.x)
}


contWt <- function(
                   q1, 
                   X, 
                   x, 
                   h,
                   profile = NULL
                   ) {         
                        # This is the function for continuous variables.
                        # q1 is the number of continuous variables.
    if (is.null(profile)) {
      profile <- X
    }
    if (q1 > 1) {
      diff.mat <- sweep(
                        X[, 1 : (q1)], 
                        2, 
                        profile[x, 1 : (q1)]
                        )
                      # This part of the function subtracts the profile of
                      # values for observation x from the profile of values 
                      # for each observation.
      kernel.wt <- kernelWeight(diff.mat / h)    
                      # This applies the kernel function to the difference 
                      # over the bandwidth, h.
      prods <- rowProds(kernel.wt)
                      # This takes the product across all q1 variables for 
                      # each observation
      prods <- ifelse(
                      is.na(prods), 
                      0, 
                      prods
                      )
                      # Sometimes rowProds has issues with zeroes in 
                      # certain columns, but it only occurs when the 
                      # product would be zero. This fixes it.
      return(prods)
    } else if (q1 == 1) {
                        # This does the same thing, but is just redone 
                        # for when there is only one continuous regressor.
      diff.mat <- X[, 1] - profile[x, 1]
      kernel.wt <- kernelWeight(diff.mat / h)
      return(kernel.wt)
    } else {
      return(1)
                        # If there are no continuous regressors, 
                        # we just multiply by one.
    }
}

ordWt <- function(
                  q1, 
                  q2, 
                  X, 
                  x, 
                  delta,
                  profile = NULL
                  ) {                                
  if (is.null(profile)) {
    profile <- X
  }  
  if (q2 > 1) {
   diff.mat <- sweep(
                     X[, (q1 + 1) : (q1 + q2)], 
                     2, 
                     profile[x, (q1 + 1) : (q1 + q2)]
                     )
                        # Columns q1 + 1 to q1 + q2 are ordered and discrete.
   kernel.wt <- delta^abs(diff.mat)
                        # The weighting function here is the bandwidth, delta,
                        # raised to the absolute value of the difference.
   prods <- rowProds(kernel.wt)
   prods <- ifelse(
                   is.na(prods), 
                   0, 
                   prods
                   )
   return(prods)
  } else if (q2 == 1) {
    diff.mat <- X[, q1 + q2] - profile[x, q1 + q2]
    kernel.wt <- delta ^ abs(diff.mat)
    return(kernel.wt)
  } else {
    return(1)
  }
}

unordWt <- function(
                    q3, 
                    X, 
                    x, 
                    lambda,
                    profile = NULL
                    ) {
  if (is.null(profile)) {
    profile <- X
  }
  if (q3 > 1) {
    kernel.wt <- lambda ^ sweep(
                                X[, (ncol(X) - q3 + 1) : ncol(X)], 
                                2, 
                                profile[x, ((ncol(profile) - q3 + 1) : ncol(profile))], 
                                FUN = '!='
                                )
                        # The last q3 columns (q1+q2+2 through q1+q2+q3+1) are
                        # unordered. The exponent will return TRUE/FALSE for 
                        # each observation, where TRUE means that the two are 
                        # not the same. If they are the same, this will be 1;
                        # if they are different, it will be lambda.
    prods <- rowProds(kernel.wt)
    prods <- ifelse(
                    is.na(prods), 
                    0, 
                    prods
                    )
    return(prods)
  } else if (q3 == 1) {
    kernel.wt <- lambda ^ (X[, ncol(X)] != profile[x, ][ncol(profile)])
                        # This compares the last regressor for each 
                        # observation to the last regressor for the 
                        # observation in question.
    return(kernel.wt)
  } else {
    return(1)
  }
}

createProfile <- function(
                          X, 
                          loc.x,
                          x.lo,
                          x.hi,
                          int         = NULL,
                          vals        = "median",
                          length.prof = 1000
                          ) {
                        # NOTE: We use int if we have an interaction between
                        # the variable of interest and another of the
                        # covariates. It takes a vector of length two,
                        # for which the first value is the location of the
                        # interacted covariate, and the second value is the
                        # location of the interaction.
  x.seq <- seq(
               from       = x.lo,
               to         = x.hi,
               length.out = length.prof
               )
  if (is.numeric(vals)) {
    X.c <- matrix(
                  rep(
                      vals, 
                      each = length.prof
                      ), 
                  nrow = length.prof
                  )
  } else {
    X.vals <- apply(
                    X,
                    2,
                    vals
                    )
    X.c <- matrix(
                  rep(
                      X.vals,
                      each = length.prof
                      ),
                  nrow = length.prof
                  )
  }
  X.c[, loc.x] <- x.seq
  if (!is.null(int)) {
    X.c[, int[2]] <- x.seq * X.c[, int[1]]
  }
  return(X.c)
}

kernelWeight <- function(
                         u,
                         kernel = "Triangular"
                         ) {
  if (kernel == "Uniform") {
    (1 / 2) * (abs(u) <= 1)
  } else if (kernel == "Epanechnikov") {
    (3 / 4) * (1 - u ^ 2) * (abs(u) <= 1)
  } else if (kernel == "Gaussian") {  
    (1 / sqrt(2 * pi)) * exp(-(1 / 2) * u ^ 2)
  } else {
    (1 - abs(u)) * (abs(u) <= 1)
  }
}

localLogitPlot <- function(
                           X.mat, 
                           Y, 
                           bw.par, 
                           nclus     = 4, 
                           q1, 
                           q2, 
                           q3, 
                           tval      = 1.96,
                           loc.x,
                           x.lo,
                           x.hi,
                           int       = NULL,
                           vals      = "median",
                           n.plot    = 1000,
                           size.line = 1.10,
                           col.line  = "blue",
                           .print    = TRUE
                           ) {
  require(ggplot2)
  k <- ncol(X.mat)
  complete <- complete.cases(cbind(Y, X.mat))
  X.mat <- X.mat[complete, ]
  Y <- Y[complete]
  profile.mat <- createProfile(
                               X.mat, 
                               loc.x,
                               x.lo,
                               x.hi,
                               int         = int,
                               vals        = vals,
                               length.prof = n.plot
                               )
  llogit.list <- localLogit(
                            X.mat,
                            Y,
                            bw.par,
                            profile.mat,
                            nclus,
                            q1,
                            q2,
                            q3,
                            tval
                            )
  change.var <- seq(
                    from       = x.lo,
                    to         = x.hi,
                    length.out = n.plot
                    )
  llogit.df <- data.frame(
                          Probability = llogit.list$p,
                          Change      = change.var,
                          Lo          = llogit.list$lo,
                          Hi          = llogit.list$hi
                          )

  p <- ggplot(
              data = llogit.df,
              aes(
                  x = Change,
                  y = Probability
                  )
              ) +
         geom_line(
                   size   = size.line,
                   colour = col.line
                   ) +
         geom_ribbon(
                     aes(
                         ymin = Lo,
                         ymax = Hi,
                   ),
               fill   = col.line,
               alpha  = 0.3,
               colour = NA
               ) +
         theme_bw()
  if(.print) {
    print(p)
  } else {
    return(llogit.plot = p)
  }
}

StandVar <- function(vec) {
  mean.vec <- mean(
                   vec,
                   na.rm = TRUE
                   )
  sd.vec <- sd(
               vec,
               na.rm = TRUE
               )
  std.vec <- (vec - mean.vec) / sd.vec
  return(std.vec)
}

###########################################################################
# Multiplot Function                                                      #
#                                                                         #
# This function is modified from the                                      #
# multiplot source code available at:                                     #
# http://www.cookbook-r.com/Graphs/Multiple_graphs_on_one_page_(ggplot2)/ #
###########################################################################

multiplot <- function(
                      ..., 
                      plotlist = NULL, 
                      file, 
                      cols     = 1, 
                      layout   = NULL
                      ) {
  require(grid)
  plots <- c(
             list(...), 
             plotlist
             )
  numPlots <- length(plots)
  if (is.null(layout)) {
    layout <- matrix(
                     seq(
                         1, 
                         cols * ceiling(numPlots/cols)
                         ),
                     ncol = cols, 
                     nrow = ceiling(numPlots/cols)
                     )
  }

 if (numPlots==1) {
   print(plots[[1]])
 } else {
   grid.newpage()
   pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
   for (i in 1:numPlots) {
     matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
    print(
          plots[[i]], 
          vp = viewport(
                        layout.pos.row = matchidx$row,
                        layout.pos.col = matchidx$col
                        )
          )
    }
  }
}

##################################
# Function to generate plot for: #
# Sender power                   #
# Target power                   #
# Number of exports              #
# Portfolio concentration (with  #
# regime type interaction)       #
#                                # 
# We will use this for multiple  #
# figures                        #
##################################

MakePlot <- function(
                     dv,
                     bw,
                     leg.pos,
                     xmat   = x.mat,
                     q1.val = q1,
                     q2.val = q2,
                     q3.val = q3
                     ) {
  plot.send <- localLogitPlot(
                              X.mat     = xmat, 
                              Y         = dv,
                              bw.par    = bw,
                              nclus     = 4,
                              q1        = q1.val,
                              q2        = q2.val,
                              q3        = q3.val,
                              tval      = 1.654,
                              loc.x     = loc.send,
                              x.lo      = min(
                                              xmat[, loc.send], 
                                              na.rm = TRUE
                                              ),
                              x.hi      = max(
                                              xmat[, loc.send], 
                                              na.rm = TRUE
                                              ),
                              vals      = vals.profile,
                              n.plot    = 1000,
                              size.line = 1.10,
                              col.line  = "blue",
                              .print    = FALSE
                              )

  plot.targ <- localLogitPlot(
                              X.mat     = xmat, 
                              Y         = dv, 
                              bw.par    = bw, 
                              nclus     = 4, 
                              q1        = q1.val, 
                              q2        = q2.val, 
                              q3        = q3.val, 
                              tval      = 1.654,
                              loc.x     = loc.targ,
                              x.lo      = min(
                                              xmat[, loc.targ], 
                                              na.rm = TRUE
                                              ),
                              x.hi      = 5,
                              vals      = vals.profile,
                              n.plot    = 1000,
                              size.line = 1.10,
                              col.line  = "blue",
                              .print    = FALSE
                              )
                          # NOTE: We can't use the maximum value of target power
                          # because we have extreme outliers that have
                          # no sufficiently close observations. However,
                          # over 99% of our observations are below 5.
  
  plot.num <- localLogitPlot(
                             X.mat     = xmat, 
                             Y         = dv, 
                             bw.par    = bw, 
                             nclus     = 4, 
                             q1        = q1.val, 
                             q2        = q2.val, 
                             q3        = q3.val, 
                             tval      = 1.654,
                             loc.x     = loc.num,
                             x.lo      = min(
                                             xmat[, loc.num], 
                                             na.rm = TRUE
                                             ),
                             x.hi      = max(
                                             xmat[, loc.num], 
                                             na.rm = TRUE
                                             ),
                             vals      = vals.profile,
                             n.plot    = 1000,
                             size.line = 1.10,
                             col.line  = "blue",
                             .print    = FALSE
                             )
  
  p.targ <- plot.targ + 
            xlab("Target Power") +
            ylab("Probability of Sanctions Success") +
            scale_x_continuous(
                               breaks = c(
                                          min(
                                              xmat[, loc.targ], 
                                              na.rm = TRUE
                                              ),
                                          5
                                          ),
                               labels = c(
                                          "Low",
                                          "High"
                                          )
                               ) 
  num.vals <- c(
                20, 
                40, 
                60, 
                80
                )
  points.vals <- (num.vals - mean(na.omit(sanctions$targ_export_variety))) /
                 sd(na.omit(sanctions$targ_export_variety))
                        # Create values for x axis.

  p.num <- plot.num + 
           xlab("Number of Exports") +
           ylab("Probability of Sanctions Success") +
           scale_x_continuous(
                              breaks = points.vals,
                              labels = num.vals
                              )
  
  p.send <- plot.send + 
            xlab("Sender Power") +
            ylab("Probability of Sanctions Success") +
            scale_x_continuous(
                               breaks = c(
                                          min(xmat[, loc.send], 
                                              na.rm = TRUE
                                              ),
                                          max(xmat[, loc.send], 
                                              na.rm = TRUE
                                              )
                                          ),
                               labels = c(
                                          "Low",
                                          "High"
                                          )
                               ) 
  
  ################################
  # Plot HH Index for democratic #
  # and autocratic targets       #
  ################################
  
  vals.dem <- vals.aut <- vals.profile
  
  vals.dem["polity"] <- 8
  
  vals.aut["polity"] <- -8
  
  plot.h.dem <- localLogitPlot(
                               X.mat     = xmat, 
                               Y         = dv, 
                               bw.par    = bw, 
                               nclus     = 4, 
                               q1        = q1.val, 
                               q2        = q2.val, 
                               q3        = q3.val, 
                               tval      = 1.654,
                               loc.x     = loc.h,
                               x.lo      = min(
                                               xmat[, loc.h], 
                                               na.rm = TRUE
                                               ),
                               x.hi      = max(
                                               xmat[, loc.h], 
                                               na.rm = TRUE
                                               ),
                               vals      = vals.dem,
                               n.plot    = 1000,
                               size.line = 1.10,
                               col.line  = "blue",
                               .print    = FALSE
                               )
  
  plot.h.aut <- localLogitPlot(
                               X.mat     = xmat, 
                               Y         = dv, 
                               bw.par    = bw, 
                               nclus     = 4, 
                               q1        = q1.val, 
                               q2        = q2.val, 
                               q3        = q3.val, 
                               tval      = 1.654,
                               loc.x     = loc.h,
                               x.lo      = min(
                                               xmat[, loc.h], 
                                               na.rm = TRUE
                                               ),
                               x.hi      = max(
                                               xmat[, loc.h], 
                                               na.rm = TRUE
                                               ),
                               vals      = vals.aut,
                               n.plot    = 1000,
                               size.line = 1.10,
                               col.line  = "blue",
                               .print    = FALSE
                               )
  
  df.demaut <- data.frame(Change      = c(
                                          plot.h.dem$data$Change,
                                          plot.h.aut$data$Change
                                          ),
                          Hi          = c(
                                          plot.h.dem$data$Hi,
                                          plot.h.aut$data$Hi
                                          ),
                          Lo          = c(
                                          plot.h.dem$data$Lo,
                                          plot.h.aut$data$Lo
                                          ),
                          Probability = c(
                                          plot.h.dem$data$Probability,
                                          plot.h.aut$data$Probability
                                          ),
                          Regime      = c(
                                          rep(
                                              "Democracy",
                                              times = nrow(plot.h.dem$data)
                                              ),
                                          rep(
                                              "Autocracy",
                                              times = nrow(plot.h.aut$data)
                                              )
                                          )
                          )
  
  p.h.demaut <- ggplot(
                       df.demaut,
                       aes(
                           x        = Change,
                           y        = Probability,
                           colour   = Regime,
                           fill     = Regime,
                           linetype = Regime
                           ) 
                       ) +
                  geom_line(size = 1.10) +
                  geom_ribbon(
                              aes(
                                  ymin = Lo,
                                  ymax = Hi,
                                  ),
                              alpha  = 0.3,
                              colour = NA
                              ) +
                  xlab("Portfolio Concentration") +
                  ylab("Probability of Sanctions Success") +
                  scale_x_continuous(
                                     breaks = c(
                                                min(
                                                    xmat[, loc.h], 
                                                    na.rm = TRUE
                                                    ),
                                                max(
                                                    xmat[, loc.h], 
                                                    na.rm = TRUE
                                                    )
                                                ),
                                     labels = c(
                                                "Low",
                                                "High"
                                                )
                                     ) +
                  scale_color_manual(
                                     values = c(
                                                "blue", 
                                                "red"
                                                ), 
                                     name   = "Regime Type"
                                     ) +
                  scale_fill_manual(
                                    values = c(
                                               "blue", 
                                               "red"
                                               ), 
                                    name   = "Regime Type"
                                    ) +
                  scale_linetype_manual(
                                        values = c(
                                                   1, 
                                                   2
                                                   ), 
                                        name   = "Regime Type"
                                        ) +
                  theme_bw()+
                  theme(legend.position = leg.pos) 
  
  #################
  # Create figure #
  #################
  multiplot(
            p.send,
            p.num,
            p.targ,
            p.h.demaut,
            cols = 2
            )
}

#########################
# Function to calculate #
# predicted values from #
# logit models          #
#########################

GetPreds <- function(
                     xmat,
                     coefs,
                     covmat,
                     var.change,
                     interaction = NULL,
                     change.lo   = NULL,
                     change.hi   = NULL,
                     change.vec  = NULL,
                     var.sq      = NULL,
                     app.fn      = "median",
                     t.diff      = NULL,
                     tval        = 1.96,
                     n.out       = 1000,
                     ci          = FALSE
                     ) {
                        # xmat is a matrix of x variables
                        # Note that xmat should include a constant.
                        # coefs is a vector of coefficient values
                        # covmat is the covariance matrix
                        # var.change is a number specifying which X
                        # is the variable of interest. This should
                        # be a number greater than 1, since a 1 would
                        # indicate the intercept.
                        # Interaction is a vector of length 2, with
                        # a first element that gives the column
                        # representing the variable with which
                        # var.change is to be interacted, and the
                        # second giving the column in which the
                        # interaction should be stored.
                        # var.sq is used if we want to include a 
                        # squared term. Just tell getPreds() which
                        # column that should be.
    k <- ncol(xmat)
    b <- coefs
    if (is.null(change.vec)) {
        change.lo <- ifelse(
                            is.null(change.lo),
                            min(
                                xmat[, var.change],
                                na.rm = TRUE
                                ),
                            change.lo
        )
        change.hi <- ifelse(
                            is.null(change.hi),
                            max(
                                xmat[, var.change],
                                na.rm = TRUE
                                ),
                            change.hi
        )
        change.vec <- seq(
                          from       = change.lo,
                          to         = change.hi,
                          length.out = n.out
        )
    }
    if (is.numeric(app.fn)) {
        xc <- app.fn
    } else {
        x.nona <- na.omit(xmat)
        xc <- apply(
                    x.nona,
                    2,
                    app.fn
        )
                        # Make sure there are no NAs when we apply the
                        # function to the matrix.
    }
    if (is.numeric(t.diff)) {
                        # t.diff is used to change the time
                        # easily, since we might care about
                        # particular values. It is a vector
                        # with two elements: the first element
                        # is the location of the original time
                        # variable, and the second is its
                        # desired value.
        t.vals <- c(
                    t.diff[2],
                    t.diff[2] ^ 2,
                    t.diff[2] ^ 3
                    )
        xc[t.diff[1] : (t.diff[1] + 2)] <- t.vals
    }
    xc.mat <- matrix(
                     rep(
                         xc,
                         each = n.out
                         ),
                     ncol = k,
                     nrow = n.out
                     )
                        # Create a matrix of constant values.
    xc.mat[, var.change] <- change.vec
                        # Vary the variable that needs varying.
    if (is.numeric(var.sq)) {
      xc.mat[, var.sq] <- change.vec ^ 2
    }
    if (length(interaction) == 2) {
      xc.mat[, interaction[2]] <- change.vec * xc[interaction[1]]
    }
    if (ncol(xc.mat) != nrow(as.matrix(b))) {
      b <- t(as.matrix(b))
    }
    xb <- xc.mat %*% b
    preds <- plogis(xb)
    if (ci) {
        deriv <- dlogis(xb)
                        # dlogis() is the derivative of plogis()
        deriv.mat <- sweep(
                           xc.mat,
                           1,
                           deriv,
                           '*'
                     )
        var.pp <- deriv.mat %*% covmat %*% t(deriv.mat)
        se.pp <- sqrt(diag(var.pp))
        ci.lo <- preds - (tval * se.pp)
        ci.lo <- ifelse(
                        ci.lo < 0,
                        0,
                        ifelse(ci.lo > 1,
                               1,
                               ci.lo)
                        )
        ci.hi <- preds + (tval * se.pp)
        ci.hi <- ifelse(
                        ci.hi < 0,
                        0,
                        ifelse(ci.hi > 1,
                               1,
                               ci.hi)
                        )
        return(
               list(
                    pred   = preds,
                    change = change.vec,
                    lo     = ci.lo,
                    hi     = ci.hi
                    )
        )
    } else {
        return(
               list(
                    pred   = preds,
                    change = change.vec
                    )
        )
    }
}

#####################################
# Function to calculate predicted   #
# probabilities, given a profile of #
# values from an x matrix           #
#####################################

LogitPlots <- function(
                       x.mat,
                       beta,
                       vcov,
                       vals.profile,
                       leg.pos,
                       .plot = TRUE
                       ) {
  p.send <- GetPreds(
                     xmat       = x.mat,
                     coefs      = beta,
                     covmat     = vcov,
                     var.change = 1,
                     app.fn     = vals.profile,
                     tval       = 1.645,
                     ci         = TRUE
                     )
  
  p.targ <- GetPreds(
                     xmat       = x.mat,
                     coefs      = beta,
                     covmat     = vcov,
                     var.change = 2,
                     app.fn     = vals.profile,
                     tval       = 1.645,
                     ci         = TRUE
                     )
  
  p.num <- GetPreds(
                    xmat       = x.mat,
                    coefs      = beta,
                    covmat     = vcov,
                    var.change = 8,
                    app.fn     = vals.profile,
                    tval       = 1.645,
                    ci         = TRUE
                    )
  
  p.h <- GetPreds(
                  xmat        = x.mat,
                  coefs       = beta,
                  covmat      = vcov,
                  var.change  = 5,
                  interaction = c(
                                  10, 
                                  9
                                  ),
                  app.fn      = vals.profile,
                  tval        = 1.645,
                  ci          = TRUE
                  )
  
  vals.aut <- vals.profile
  vals.aut[10] <- -8
                          # Construct a profile for which Polity is
                          # set to -8.
  
  p.aut <- GetPreds(xmat        = x.mat,
                    coefs       = beta,
                    covmat      = vcov,
                    var.change  = 5,
                    interaction = c(
                                    10, 
                                    9
                                    ),
                    app.fn      = vals.aut,
                    tval        = 1.645,
                    ci          = TRUE
                    )
  
  ###################################
  # Create data frames for plotting #
  ###################################
  
  df.send <- data.frame(
                        Change      = p.send$change,
                        Hi          = p.send$hi,
                        Lo          = p.send$lo,
                        Probability = p.send$pred
                        )
  
  df.targ <- data.frame(
                        Change      = p.targ$change,
                        Hi          = p.targ$hi,
                        Lo          = p.targ$lo,
                        Probability = p.targ$pred
                        )
  
  df.num <- data.frame(
                       Change      = p.num$change,
                       Hi          = p.num$hi,
                       Lo          = p.num$lo,
                       Probability = p.num$pred
                       )
  
  df.h <- data.frame(
                     Change      = c(
                                     p.h$change,
                                     p.aut$change
                                     ),
                     Hi          = c(
                                     p.h$hi,
                                     p.aut$hi
                                     ),
                     Lo          = c(
                                     p.h$lo,
                                     p.aut$lo
                                     ),
                     Probability = c(
                                     p.h$pred,
                                     p.aut$pred
                                     ),
                     Regime      = c(
                                     rep(
                                         "Democracy",
                                         times = length(p.h$change)
                                         ),
                                     rep(
                                         "Autocracy",
                                         times = length(p.aut$change)
                                         )
                                     )
                     )
  
  ##################
  # Generate plots #
  ##################
  
  p.send <- ggplot(
                   df.send,
                   aes(
                       x = Change,
                       y = Probability
                       ) 
                   ) +
                   geom_line(
                             size   = 1.10,
                             colour = "blue"
                             ) +
                   geom_ribbon(
                               aes(
                                   ymin = Lo,
                                   ymax = Hi,
                                   ),
                               alpha  = 0.3,
                               fill   = "blue",
                               colour = NA
                               ) +
                   xlab("Sender Power") +
                   ylab("Probability of Sanctions Success") +
                   scale_x_continuous(
                                      breaks = c(
                                                 min(
                                                     x.mat[, 1], 
                                                     na.rm = TRUE
                                                     ),
                                                 max(
                                                     x.mat[, 1], 
                                                     na.rm = TRUE
                                                     )
                                                 ),
                                      labels = c(
                                                 "Low",
                                                 "High"
                                                 )
                                      ) +
                   ylim(
                       0,
                       1
                       ) +
                   theme_bw()
  
  p.targ <- ggplot(
                   df.targ,
                   aes(
                       x = Change,
                       y = Probability
                       ) 
                   ) +
                   geom_line(
                             size   = 1.10,
                             colour = "blue"
                             ) +
                   geom_ribbon(
                               aes(
                                   ymin = Lo,
                                   ymax = Hi,
                                   ),
                               alpha  = 0.3,
                               fill   = "blue",
                               colour = NA
                               ) +
                   xlab("Target Power") +
                   ylab("Probability of Sanctions Success") +
                   scale_x_continuous(breaks = c(
                                                 min(
                                                     x.mat[, 2], 
                                                     na.rm = TRUE
                                                     ),
                                                 max(
                                                     x.mat[, 2], 
                                                     na.rm = TRUE
                                                     )
                                                 ),
                                      labels = c(
                                                 "Low",
                                                 "High"
                                                 )
                                      ) +
                   ylim(
                       0,
                       1
                       ) +
                   theme_bw()
  
  p.num <- ggplot(
                  df.num,
                  aes(
                      x = Change,
                      y = Probability
                      ) 
                  ) +
                  geom_line(
                            size   = 1.10,
                            colour = "blue"
                            ) +
                  geom_ribbon(
                              aes(
                                  ymin = Lo,
                                  ymax = Hi,
                                  ),
                              alpha  = 0.3,
                              fill   = "blue",
                              colour = NA
                              ) +
                  xlab("Number of Exports") +
                  ylab("Probability of Sanctions Success") +
                  scale_x_continuous(
                                     breaks = c(
                                                min(
                                                    x.mat[, 8], 
                                                    na.rm = TRUE
                                                    ),
                                                max(
                                                    x.mat[, 8], 
                                                    na.rm = TRUE
                                                    )
                                                ),
                                     labels = c(
                                                "Low",
                                                "High"
                                                )
                                     ) +
                  ylim(
                       0,
                       1
                       ) +
                  theme_bw()
  
  p.h <- ggplot(
                df.h,
                aes(
                    x        = Change,
                    y        = Probability,
                    colour   = Regime,
                    fill     = Regime,
                    linetype = Regime
                    ) 
                ) +
                geom_line(size   = 1.10) +
                geom_ribbon(
                            aes(
                                ymin = Lo,
                                ymax = Hi,
                                ),
                            alpha  = 0.3,
                            colour = NA
                            ) +
                xlab("Portfolio Concentration") +
                ylab("Probability of Sanctions Success") +
                  scale_x_continuous(breaks = c(
                                                min(
                                                    x.mat[, 5], 
                                                    na.rm = TRUE
                                                    ),
                                                max(
                                                    x.mat[, 5], 
                                                    na.rm = TRUE
                                                    )
                                                ),
                                     labels = c(
                                                "Low",
                                                "High"
                                                )
                                     ) +
                  scale_color_manual(
                                     values = c(
                                                "blue", 
                                                "red"
                                                ), 
                                     name   = "Regime Type"
                                     ) +
                  scale_fill_manual(
                                    values = c(
                                               "blue", 
                                               "red"
                                               ), 
                                    name   = "Regime Type"
                                    ) +
                  scale_linetype_manual(
                                        values = c(
                                                   1, 
                                                   2
                                                   ), 
                                        name = "Regime Type"
                                        ) +
                  ylim(
                       0,
                       1
                       ) +
                  theme_bw()+
                  theme(legend.position = leg.pos) 
  if (.plot == TRUE) {
    multiplot(
              p.send,
              p.num,
              p.targ,
              p.h,
              cols = 2
              )
  } else {
    return(
           list(
                send = p.send,
                num  = p.num,
                targ = p.targ,
                h    = p.h
                )
           )
  }
} 
  
CorrPred <- function(
                     pred,
                     actual,
                     thresholds = seq(
                                      0.1,
                                      0.9,
                                      by = 0.1
                                      )
                     ) {
  correct.df <- data.frame(
                           thresholds = thresholds,
                           prop.corr  = NA,
                           false.pos  = NA
                           )
  for (i in 1 : length(thresholds)) {
    pred.val <- ifelse(
                       pred >= thresholds[i],
                       1,
                       0
                       )
    correct.df$prop.corr[i] <- sum(pred.val == actual) / length(actual)
    correct.df$false.pos[i] <- sum(pred.val > actual) / length(actual)
  }
  return(correct.df)
}

CorrPredPlot <- function(
                         y,
                         x1,
                         x2,
                         b1,
                         b2
                         ) {
  pred1 <- na.omit(
                   cbind(
                         y,
                         plogis(x1 %*% b1)
                         )
                   )
  pred2 <- na.omit(
                   cbind(
                         y,
                         plogis(x2 %*% b2)
                         )
                   )
  corr1 <- CorrPred(
                    pred1[, 2], 
                    pred1[, 1]
                    )
  
  corr2 <- CorrPred(
                    pred2[, 2], 
                    pred2[, 1]
                      )
  table.corr <- data.frame(
                           Threshold  = corr1[, 1],
                           Prop.Col1  = corr1[, 2],
                           Prop.Col5  = corr2[, 2]
                           )
  corr.df <- melt(
                  table.corr,
                  id.vars = "Threshold"
                  )
  p <- ggplot(
              corr.df,
              aes(
                  x        = Threshold,
                  y        = value,
                  colour   = variable,
                  linetype = variable
                  )
              ) +
         geom_line(
                   size  = 1.25,
                   alpha = 0.8
                   ) +
         ylab("Proportion Correctly Predicted") +
         scale_colour_discrete(
                               name   = "Model",
                               labels = c(
                                          "Column 1",
                                          "Column 5"
                                          )
                               ) +
         scale_linetype_discrete(
                                 name   = "Model",
                                 labels = c(
                                            "Column 1",
                                            "Column 5"
                                            )
                                 ) +
         theme_bw()  
  return(p)
}
